One Parameter Scaling Theory for Stationary States of Disordered Nonlinear Systems 
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We show, using detailed numerical analysis and theoretical arguments, that the normalized par- 
ticipation number of the stationary solutions of disordered nonlinear lattices obeys a one-parameter 
scaling law. Our approach opens a new way to investigate the interplay of Anderson localization 
and nonlinearity based on the powerful ideas of scaling theory. 
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Introduction: Wave propagation in naturally occurring 
or engineered complex media is an interdisciplinary field 
of research that addresses systems as diverse as classical, 
quantum and atomic-matter waves. Despite this diver- 
sity, the wave nature of these systems provides a common 
framework for understanding their transport properties 
and often leads to new applications. One such charac- 
teristic is wave interference phenomena. Their existence 
results in a complete halt of wave propagation in random 
media, which can be achieved by increasing the random- 
ness of the medium. This phenomenon was predicted 
fifty years ago in the framework of quantum (electronic) 
waves by Anderson [1 and its existence has been con- 
firmed in recent years in experiments with classical [21410] 
and matter waves [TTJ [T^] . 

In many of these experiments, the appearance of non- 
linearities, induced either due to the nonlinear Kerr-effect 
(in the framework of nonlinear wave propagation in dis- 
ordered photonic lattices) [8HTU] or due to atom-atom in- 
teractions (in atomic transport of Bose-Einstein Conden- 
sates in optical lattices) [TTJ [T^] , might affect drastically 
the phenomenon of Anderson localization. An important 
question is therefore how the interplay between disorder 
and nonlinearity might complement, frustrate, or rein- 
force each other [T3HI6] . This notion is not only of exper- 
imental importance, but it raises a number of unsolved 
theoretical questions as well. The theoretical study of 
localization in random nonlinear lattices has been ad- 
vanced using several approaches; including the studies of 
transmission [17], wavepacket dynamics |18j . and station- 
ary solutions [TM2"T] . 

In this Letter, we approach the interplay of nonlin- 
earity with disorder from a different perspective, namely 
we develop a scaling theory for localization phenomena 
in nonlinear random media described by the Discrete 
Nonlinear Schrodinger Equation (DNLSE) . Scaling ideas 
played a major role in understanding various proper- 
ties of linear disordered systems, including the structure 
of their eigenstates However, solutions of the 

DNLSE, for sufficiently strong nonlinearity, have noth- 
ing in common with the solutions of the linear problem. 
Indeed, the number of solutions of the DNLSE is gener- 



ally much larger than the number of eigenstates of the 
corresponding linear problem and, in particular, many 
solutions appear outside the spectrum of the linear sys- 
tem. It is therefore quite remarkable that the scaling 
ideas can be extended to the nonlinear case. Specifically, 
we find that the rescaled participation number p^ (x) of 
the stationary solutions of the DNLSE of lattice size N 
and nonlinearity Xi obeys a one-parameter scaling, i.e. 
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Above /3 is a universal function of p^ alone, which is 
independent of any microscopic properties of the system 
under investigation, and (• • •) denotes an averaging over 
disorder realizations and over states within a small fre- 
quency window. The participation number £jv is |26j 
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and is proportional to the effective number of nonzero 
components i/j n of a stationary solution of the DNLSE. 
£jy is the participation number for some reference ensem- 
ble, chosen such that it supports the most extended states 
for a specific lattice topology and nonlinearity. It will be 
argued below that t^f is proportional to the system size 
N. Note, however, that it would not be accurate to re- 
place £™ by N because the coefficient of proportionality 
between the two lengths depends weakly on x- Eq.Q is 
confirmed in the following via detailed numerical simu- 
lations with quasi one-dimensional (ID) disordered sys- 
tems described by Banded Random Matrix (BRM) mod- 
els and is supported by theoretical arguments [27] . 

Mathematical Model: We consider a class of random 
systems described by the time-independent DNLSE: 



^ H nm ip m + x\i>n\ 2 ipn = uvb r , 
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where u) are the frequencies of the stationary solutions 
and tp n is their amplitude at site n. The connectivity 
matrix H defines the topology of the sample. In the case 
of strictly ID disordered systems, it is a three-diagonal 
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matrix with H n>n ±i = — 1 while H nn are random in- 
dependent variables given by some distribution. In our 
simulations below, we consider the challenging case where 
H belongs to a BRM ensemble, having in mind quasi-l-D 
systems with b propagating modes [531 HI] • 

The BRM ensemble is defined as a set of real symmetric 
N x N matrices with elements H nm = for \n — m\ > 6, 
while inside the bandwidth b they are independent ran- 
dom variables given by a Gaussian distribution of mean 
zero and fixed variance 123 



(Hnm) — 0, 



b(2N- 6 + 1) 



With this normalization, the eigenvalues of H (in the 
limit of large N, b) are located in the interval (—2, 2) [23] . 

Numerical Method: The stationary solutions of the 
nonlinear Eq.([3]) were found numerically by utilizing a 
continuation method approach. Starting with the linear 
modes of H as an initial guess, we take a small step 
in x (<^X ~ 10 4 ), and solve the nonlinear system of 
Eq.([3]). This is achieved by minimizing a multivariable 
(N + l)-dimensional vector function F n = J2 m H n m4>m + 
X\i>n\ 2 ipn-ui} n ,n= 1,2, ..JV and F N +i = J2 n l^n| 2 -l- 
Using this method we find ({ip n };cj) with tolerance 10~ 8 . 
The resulting solution then becomes the initial guess for 
the nonlinear solver for the next step in x- 

An "evolution" of a representative Anderson localized 
mode as nonlinearity x increases is reported in Fig. [T] 
We observe that while initially (small x values) its shape 
remains unaffected, eventually it starts delocalize until it 
spreads over the whole sample. The degree of derealiza- 
tion as a function of x 1S reflected in the behavior of the 
participation number 6v(x) (lower panel of Fig. [I]). We 
want to investigate how the average participation num- 
ber (£/v) of the stationary modes of Eq.([3| is affected by 
nonlinearity. The averaging (•••), has been performed 
over solutions with corresponding uj being inside a small 
frequency interval (below w € [—1,1]) such that the na- 
ture of wavefunctions is statistically the same. For better 
statistical processing, a number of disorder realizations 
has been used, such that the total number of obtained 
stationary solutions is at least 10 4 . 

Nonlinear Localization Length: We start our analysis 
by introducing the asymptotic localization length \ x de- 
fined via the participation number £/v : 



lim (6v(x))- 

N— >oo 
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In Fig. [2^l we report some representative data for the 
participation number (£n(x))j as a function of the system 
size N. From these plots we extract the saturation value 
A x . The resulting data are summarized in Fig. [2Jd by 
referring to the rescaled localization length A x /Ao, where 
Aq is the localization length given by Eq. ^ for the linear 
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FIG. 1: Upper panel: Parametric evolution as the nonlinear- 
ity x increases of a stationary solution of the BRM, Eq.|3|, 
with 6 = 4. For \ — the mode is exponentially localized 
while as x increases it becomes delocalized over the whole 
lattice. Lower panel: The participation number £jv(x) (solid 
black line) as a function of X- The coloring reflects the wave- 
function intensity shown in the upper panel. 



(i.e. x = 0) system. We find that 
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1 for x < X* 
v/X for X» X* 
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where x* ~ 0.3 [18] . A simple interpolation formula 
which agrees with our numerical results (see Fig.^p) is: 



A x = AqVTT 



aoX 



(7) 



where the fitting parameter a was found to be a ~ 3. 
While interpreting the above equations, it is crucial to 
keep in mind that all the lengths depend on frequency ui 
- when comparing these lengths for different values of the 
nonlinearity Xi one should keep u> (approximately) fixed. 

The following heuristic argument provides some un- 
derstanding of the behavior of A x in the two limiting 
cases. First, let us note that the inverse of the partici- 
pation number ^ (%) is proportional to the interaction 
energy stored in the system described by the DNLSE, i.e. 
Emt = f J2 n IV'nl 4 • For x = the localization length A x 
is equal to Aq, while due to normalization the correspond- 



ing stationary solutions of Eq. (131) are ip\ 
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When x increases, Ei nt grows within the localization 
length Ao, but oj is assumed to be approximately fixed. 
This increase in E- mt is compensated by further spread- 
ing of the wavefunction beyond Ao- The first question 
is: what is the value of the nonlinearity strength \* , for 
which the spreading beyond Ao becomes significant? An 
estimation is achieved by comparing the interaction en- 

ergy E int = (x/2)EJ^ 0) | 4 « (x/2)A (l/A 2 ) = X /2X 
stored in the "localization box" of size Ao to the cor- 
responding mean level spacing Aa ~ 1/Ao. From this 
comparison, we get x* ~ 1- 
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FIG. 2: (a) The participation number £jv(x) vs - the system 
size N, for various nonlinearity strengths: \ = (°),X = 
0.1 (A),x = 1 (0),X = 10 ( V ),X = 50 (D),x = 100 (>). The 
asymptotic localization length A x is evaluated by a direct fit 
of the saturation value of fiv(x) hi the limit iV — > oo; (b) The 
extracted A x versus the nonlinearity strength x- The black 
dashed line has slope 0.5 and is drawn in order to guide the 
eye. The red dashed line is the best fitting curve Eq.Q. 

Next, we consider the limit x ^ X* i when A x 3> Ao. 
In this case, the wavefunctions ip 1 ^ ~ lyf\ spread 
over A x /Ao "localization boxes". We make the following 
self-consistent argument: The total interaction energy 
stored is E int = ( X /2) £ n |^ x) | 4 « (x/2)A x (l/A 2 ) = 
x/2A x . The interaction energy per box is therefore 
(x/2A x )(Ao/A x ) = (x/2)Ao/A x . However, from the pre- 
vious considerations we know that one localization box 
can "resist" an energy A\ ~ 1/Ao- A self-consistency 
condition gives A x ~ Ao_>/x, which agrees with the 
X ^ X* asymptotic in Eq. (pi) . 

Reference Ensemble: The most ergodic stationary so- 
lutions of Eq.([3| correspond to connectivity matrices H, 
taken from the GOE [55]. These solutions will spread 
over the entire system, of size N. However, this does not 
mean (fif = A, since ^4 might have some (oscillatory) 
structure. We therefore assume that 

etf(x) = Na( X ) where 1/3 < a( X ) < 1. (8) 



The lower value of a(x) is achieved in the limit of X —> 0, 
where (£|f) = A/3 [28] . The fact that (gjf) is 3 
times less than the system size is due to Gaussian fluc- 
tuations in the components tp n . The other limiting case 
of x ~^ 00 corresponds to a = 1. Indeed, strong non- 
linearity favors a completely uniform iV'nl 2 (again, for 
fixed and moderate oj). The following argument, simi- 
lar in spirit to the "maximum entropy" ansatz, can pro- 
vide some understanding. Let us denote the compo- 
nents of a stationary solution by a random vari- 
able x. Assume that the variable x follows a distribution 
V(x), "as random as possible" but with the constraint 
Af = J dxx 2 V{x) = 1/N dictated by normalization of the 
wavefunction. For x = 0, it was shown that the "most 
random" distribution is the Gaussian. It can be found by 
maximizing the entropy S = — J dxV(x) InV(x) with the 
above constraint [29]. For x ^ 0, one however has also to 
consider the increase of interaction energy E- m t. The en- 
tropy favors a distribution V(x), "as random as possible" 
for the values of x at various sites; e.g. the distribution 
of the "weight" when all of the probability is on a single 
site is as likely as the configuration with equal values of 
all the sites. The energy however favors a uniform distri- 
bution - clearly, a "single site weight" will require huge 
energy. Thus, the correct quantity to minimize is the 
"free energy" functional F[V] = -S + E int +/3j\f. We get 



V(x) = Cexp(-/3a; 2 - (l/2)xa; 4 ) 
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where C and /? are defined from the constraints j\f = 1 /N 
and / V(x)dx = 1. For x = the Gaussian distribution 
is recovered, whereas for x - > oo the distribution turns 
(5-like around x = 1/ y/N. 

One Parameter Scaling Ansatz: We are now equipped 
to formulate a scaling theory for the stationary solutions 
of Eq.([3|. The scaling ansatz of Eq.Q is equivalent to 
postulating the existence of a function f(x) such that 
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In the delocalized limit x 3> 1, the function f(x) ap- 
proaches the saturation value 1. On the other hand, 
when N — > oo (i.e. x <C 1) we have that (6v(x)) ~ * ^x> 
thus f(x) — x. We have numerically tested Eq.([To|, for 
the DNLSE using a BRM ensemble for the connectiv- 
ity matrix H. Various values of N and b in the ranges 
50 < N < 800 and 1 < b < N/2 were used in the analy- 
sis. The numerical data are reported in Fig.[3j confirming 
nicely the scaling ansatz of Eq.(10). Finally, it is reassur- 
ing that Eq. ( 10 ) in the limit x = 0, recovers the scaling 
relation found for linear disordered lattices [2"3Tf2"5"] . 

It is clear that "strong" Discrete Breathers (DB) lo- 
calized at few sites are excluded from our above consid- 
erations. In fact, such solutions correspond to large fre- 
quency values lo ~ Xi which for very large x are always 
outside any fixed small frequency window over which our 
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FIG. 3: Main panel: Participation ratio £jv(x)/£5v f (x) f° r the 
BRM model vs. the rescaled parameter A x /^ f (x) f° r various 
(TV, 6)-values (for each TV we have used at least 12 different 
b values in the interval 1 < b < N/2). Colors correspond 
to nonlinear values of: x = (black), x = 0.1 (red), x = 1 
(orange), \ — 10 (blue), and x — 100 (green). The solid 
line is the best fit with the function y = 0.98x/(l + 0.98a;). 
Upper inset: The same data as in the main panel, reported 
using the scaling variables {£n)/N and Xq/N. Lower Inset: 
The numerically evaluated a(x) vs. \ f° r the nonlinear GOE 
reference ensemble. The size of the matrix H is N = 800. In 
the limit N — > oo and for \ 3> 1, the function a(x) — > L 



scaling analysis is performed. The analysis of DB can 
be done separately, since they originate from the clean 
system with large x an d w, in which disorder has only 
a minor effect. Thus, all lengths appearing in Eq.(10) 



are (by their definition) approximately the same, i.e. 
Cw f ~ 6v(x) ~ <\jcj 011 t ne or der of the lattice spac- 
ing. This is a rather trivial limit, where Eq.(10) is again 
expected to be approximately valid. 

Conclusions: We presented a one-parameter scaling 
theory for the average participation number of the sta- 
tionary solutions of low-dimensional disordered nonlinear 
systems described by the DNLSE. Via numerical analy- 
sis and theoretical considerations, we have established 
Eq.([T]) which allows us to conclude that changing disor- 
der strength, nonlinearity, and frequency, in the way de- 



scribed by Eq.(10), would not change the (average) spa- 



tial extent of the stationary states of the DNLSE. The 
one-parameter scaling theory presented here is a power- 
ful approach in the quest of understanding the interplay 
between nonlinearity and disorder. Although the focus 
of this Letter was on the structure of stationary modes, 
our approach can further be used to understand various 
other observables [3D]. 
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